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Abstract. - We present exact solutions of the non-linear BCRE model for granular avalanches 
without diffusion. We assume a generic sandpile profile consisting of two regions of constant but 
different slope. Our solution is constructed in terms of characteristic curves from which several 
novel predictions for experiments on avalanches are deduced: Analytical results are given for 
the shock condition, shock coordinates, universal quantities at the shock, slope relaxation at 
large times, velocities of the active region and of the sandpile profile. 



Introduction and Model. - The study of avalanches and surface flows in granular materials 
has attracted much attention recently, both from a theoretical Q and an experimental point 
of view [Q. A simple model, thought to capture some of the essential phenomena, has been 
proposed in || |5j. It is based on the assumption that a strict separation between rolling 
grains and static grains can be made. Coupled dynamical equations for these two species, 
based on phenomenological arguments, can then be written. Calling R the local density of 
rolling grains and h the height of static grains, the simplest form of the BCRE equations read: 

H t = ~jRH x , (1) 
R t = R x +jRH x , (2) 

where H is the height of static grains, counted from the repose slope of angle 8 r : h(t,x) = 
H(t,x) + xtan(0 r ) (the heap is sloping upwards from left to right). In the above equations, 
the units of lengths and time are chosen such that the (downhill) velocity of grains is v = 1 , 
while H and R are counted in units of the grains diameter. The term "/RH X describes the 
conversion of static grains into rolling grains if H x > 0, or vice versa if H x < 0. 7 is a grain 
collision frequency, typically of the order of 100 Hz. 

Many important phenomenon are left out from the above description, and can be included 
by adding more terms. For example, diffusion terms (such as D\R XX or E>iH xx , describing, 
e.g. non local dislodgement effects) will generically be present, and qualitatively change the 
structure of the solutions |5j . Another aspect not described by the linear form of the conversion 
term above is the expected saturation of rolling grains with time, rather than the exponential 
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growth predicted by Eq. (||) for a constant positive slope H x . Non linear saturation terms, as 
well as a dependence of the velocity of the rolling grains on R, are thus expected in general, 
and can lead to important differences with the above equations || 

Recently, these equations has been studied by Mahadevan and Pomeau (mp) ||. They 
found a conservation law, which relates the solutions R(t, x) and H(t, x) in a frame moving 
with the velocity of the grains. From this law, they concluded that the BCRE equations have 
characteristics that are straight lines, along which both R(t,x) and H(t,x) are constant. 
Independently of the initial profile Hq(x), they found that a shock forms at time t s = 
— 1/ ('yR'o max ) with i?g max is the maximum (in absolute value) of the initial gradient of rolling 
grains. Whereas our exact solution fulfills the same conservation law, our results for the 
characteristics and the shock time disagree with the results of MP. As we will discuss below, 
the reason for this disagreement is their implicit assumption of a very restrictive relation 
between the initial profiles Ro(x) and Hq(x). 

Characteristic coordinates. - The general basis of the method Q we used to solve Eqs. 
(0J|) consists in a replacement of the original equations by an equivalent system of four partial 
differential equations for the functions t, x, R and H, but now considered as functions of new 



coordinates /i and v, which will be defined below. ( 1 ) These new equations will be particularly 
simple inasmuch as each equation has derivatives with respect to either \i or v, though the 
mapping between the coordinates (t,x) and (fJ,,v) will be in general complicated. To define 
the characteristic coordinates (fJ.,v), we have to specify first the characteristic curves of the 
system (0,||). For practical reasons, we introduce new functions u(t,x) = 1 — R(t,x)/a and 
v(t, x) = (a + x — x) — H(t, x))/a instead of H(t, x) and R(t, x). For this new functions 
the differential expressions become 

Li[u,v] = -u t — ^a{l-u)u x + v t +^a{l-u)v x -^{l-u)=Q, (3) 
L 2 [u,v] = u t + [-1 +7«(1 - u))u x -70(1 - u)v x + 7(1 - u) = 0. (4) 

Both operators L± and L2 contain linear combinations of the type aut+bu x of the derivatives of 
u (and the same holds for v). This combination means that u is differentiated in the direction 
given by the ratio t/x = a/b. Since the coefficients a and b differ for u and v and also for L\ 
and L2, the functions u and v arc differentiated in each of the operators in different directions 
in the (t, x) plane. Notice that the directions depend also on u itself, and therefore on the 
solution under consideration, which is a typical feature of non-linear systems. As noted above, 
our goal is to find equivalent differential equations of which each contains derivatives in only 
one (local) direction corresponding to one of the new coordinates u and v. Therefore we take 
a linear combination L — \\Li + X2L2 of the operators in Eqs. ([§||) such that the derivatives 
of u and v in L combine to derivatives in the same direction, which is called a characteristic 
direction. Moreover we assume that these local directions change smoothly as functions of t 
and a;, and are given by the tangential vectors (^(tr), x a (a)) of a smooth path (t(a), x(a)) 
with a as parameter. Considering the functions u and v along this path, they depend only on 
a and we have, e.g., u a — u t t a + u x x a . Using these conditions, we obtain four homogeneous 
linear equations for the coefficients Ai and A2 with coefficients depending on t, x, u, v and their 
derivatives with respect to a. For non-trivial solutions all possible determinants of the matrix 
of these coefficients have to vanish, leading to three independent equations or characteristic 
relations (cr). The first one can be written as a quadratic equation for the local direction 
C = x a /t a of differentiation, the solution of which are: C+ = ~ 1 an d C- = jcx(l— u). Now, for a 

( x )The theory used here is actually more general and can be used in the presence of non-linear 
saturation terms or for ripple models [pi. 
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fixed solution u, the equations dx/dt — £ + and dx/dt — C- are ordinary differential equations, 
which define two families of paths with the starting position xq at t = as parameter. These 
families of paths are the characteristics C+ and C_ of the system From a physical point 

of view, they are simply the paths along which R(x, t) (C+) and H{t, x) (£_) evolves with time. 

The new curved coordinate frame (/z, v) is now defined such that the two one-parametric 
families of characteristics are mapped by the coordinate transformation on an usual Cartesian 
coordinate frame in the (/x, t/)-plane, i.e., along the characteristics the coordinate functions 
n(t,x) and v(t,x), respectively, are constant. Here we have chosen to map the line t = on 
the line given by /i = —v. In terms of the new coordinates we find 

x lJ +t l/ = 0, Xfj, — ja(l — u)tfji, = 0. (5) 

Now we make use of another CR, which evaluated along C+ and C_ by identifying a with v 
and /j,, respectively, yields the conditions 

u v + 7a(l — u)v u + 7(1 - u)t v = 0, it M - 1^ + 7(1 - ii)t M = 0. (6) 

These equations together with the Eqs. (||) form the desired set of four equations mentioned 
before. Every solution of this new system satisfies the original Eqs. since the Jacobian 

tvXp — t^Xu ~ 1 + jR(fj,, v) of the coordinate map does not vanish due to 7-R(/i, v) > 0. 

General solution. - Before we can construct a solution to the equivalent system (^,^|), we 
have to specify initial data along the line /i = —v corresponding to t = 0. We choose an 
general profile Ho(x), perturbed at t = by a uniform 'rain' of rolling grains: Rq(x) = a. In 
terms of the new coordinates, the initial conditions become to(^) = 0, Xo(fi) = —fi, uo(fi) = 0, 
vo(ii) = — (/j, + Ho(—fi))/a. By introducing the function A(fi,v) = — 1 — 7a(l — u(fi,v)) 1 
one can show that the problem of solving the system given by Eqs. (||,||) can be reduced to 
the task of finding a solution to the equation A y = ^H' Q {y){l + 1/A), with initial condition 
A(/i, — fi) = —1 — ja. The solution of this equation can be simply expressed in terms of the 
so-called Lambert function W jll| : 

A(fi,u) = -l-W{ajexp[a~/ + 7 (H (-ti)-H {v))]}. (7) 

With this solution at hand, the solution to the system is determined by 

, s r ds 

HM> v) = \ -it, 7 

J-u A(s,i/) 



where we have expressed already the original fields v) and H(fi, v) in terms of the functions 
u and v. To get the fields as functions of t and x, one has to invert the coordinate map. This 
can be done by using x) = — t — x and integrating the equation for t(jj,, v) to obtain also 
y(t,x). As announced before, the height profile H(t,x) = Ho(i/(t, x)) turns out to be constant 
along the characteristics C_ . 

Generic shape for H(t, x). - In the following we will consider a situation which is generic for 
sandpile surfaces. Suppose that one starts with a sandpile profile, which consists of two regions 
with constant but different slopes matching with a kink at x — 0, and again with a constant 
amount of rolling grains. The slopes may be either larger or smaller than the angle of repose 
8 r . If we denote the slope to the right (left) by 8 r + > (9 r +9<), we have Hq(x) = > x for 
x > and H (x) = 9 K x for x < 0. In the case of a piecewise constant H' (x) one can integrate 



= —fl — v 



" A M (s,^)d s 



-« 7#o(-*) 
1 + A(/i, v) 

7 



x(p,v) = -n-t(n,v) 
H(n,u) = H (p), (8) 
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the equation for t(pt, v) easily as can be seen from Eq. (@). The structure of Eq. (Q) suggests 
to distinguish between three regions given by fi > 0,f < (I), /j,, v < (II) and /.i < 0, v > 
(III). |Q] In regions I and III one can find the explicit expression i/(t, x) = x + j(l — e l8t ) with 
9 = 9 < (I) or 9 = 6 > (III), i.e., the characteristics C_ are in these regions simple exponential 
curves. As a consequence, no shocks can appear in these two regions and the corresponding 
solutions are particularly simple: 

Ri(in)(t,x) = ae< e <l» t , H I(ni) (t,x) = H Q (x) + a-R I(III) (t,x). (9) 

The boundaries of the regions I and III in real space (t, x) are given by the conditions x < — t 
and x > x\(t) — j-(e l9>t — 1) corresponding to the [i = and v = characteristics, 
respectively, see Fig. |l|. The boundary for region I has an obvious physical meaning: The 
information that there is a kink at x — can only propagate to the left with the velocity of 
the moving grains, which is 1 in our rescaled units. Moreover, it is important to note that the 
'uphill' velocity with which the kink moves is only equal to 7a at small times, before growing 
exponentially. As discussed in the introduction, this growth eventually saturates, as does the 
value R, or else the characteristic C_ quickly reaches the edge of the pile. 

The range of x in between the above two regions corresponds to the intermediate region II. 
Within this range one can obtain only an implicit solution for the coordinate map v(t,x). It 
reads 

A(-a; - t, u) - A(0, u) A(0, v) + 1 + cry! . , 

where A(/x, v) = — 1 — W {cry exp[«7 — "/(9 > fi + 9<v)]} as follows from Eq. (uh. The shape 
of R(t, x) and H(t, x) can be obtained directly from the last two equations of (g|). In general, 
Eq. ( [Toj ) has to be solved numerically although several results can be obtained in an analytic 
way. It turns out that the solutions of Eq. ( |l0| ) fall into two qualitatively different classes, 
according to the values of (3 = 9 > /9 < and 6 < : for j3 > 1 — 07 or 6 < < 0, both R(t,x) and 
H(t, x) remain continuous for all times, while for (3 < 1 — a-f and 9 < > 0, the solutions develop 
a discontinuity in R(t,x) and H(t,x) beyond a finite shock time t s . This must be contrasted 
with MP, since in the present case Ro(x) = a, they predict that shocks are absent for all times. 

Examples. - The characteristics resulting from numerical solutions of Eq. ( [To| ) have been 
plotted in Fig. |l|. The left part of this figure has been obtained for 9 > > and 9 < < 0, 
corresponding to (3 < 0. In this case, the characteristics are more and more 'diluted' as time 
increases, and therefore never cross - no shock. In the limit of large times, the argument of 
the Lambert W function becomes very large. Using the first two terms of the asymptotic 
expansion of W we get i/(t, x) = [—{it + ln(x + t^j)/(70<)]/(/9 - !)■ The corresponding 
expression for R(t,x) and H(t,x) can be obtained from Eq. (||). A particularly interesting 
quantity to look at is the local slope at, say, x = 0. In this limit the slope is negative and 
decays with time as H x (t,x — 0) — l/(7/3t).|( 5 )| It means that the 'true' slope h x actually 
relaxes to the angle of repose 9 r for very large time. If L is the size of experimental system, 
then C- reaches the boundary of the system at a time t* such that L « j-e lB>t . One should 
therefore measure a final slope h x m 9 r + 9 < j \i\(9 > L / a) smaller than the repose angle. This 
result is consistent with the qualitative discussion of Boutreux and de Gennes for a similar 
situation | |l3fl . 

( 2 ) The region where /1, v > turns out to be mapped on the half space with t < and is therefore 
not of physical interest. 

( 3 ) Note that this t~ x relaxation of the slope has also been obtained in jl2j within a very different 
model. 
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Another experimentally important quantity is the velocity Vr of the "active" region. Fol- 
lowing ]5|], this region can be defined by the condition R(t,x) > i? m in, where R m i n is a small 
threshold, vr is then given by the slope of the curves of constant R(t,x), which tends to 
a constant in the large t limit as can be seen in Fig. |l](a). The asymptotic analysis yields 
vr = /3/(l — /3). Since /3 < 0, —1 < vr < 0, and the avalanche proceeds downhill, but slower 
than the grains themselves. This is an effect of the non-linear term in the BCRE equations 
since the linearized theory yields vr = — 1 ||. 




Fig. 1. - Characteristics C_ for the cases (a) 0> = 0.1, 0< = -0.1 and (b) 6> = -0.1, 0< = 0.1. 
For both cases we have taken 7 = 0.5, a = 0.4. The dashed lines represent the boundaries between 
the different regions I, II and III explained in the text. The bold line in (b) is the envelop of the 
characteristics. It presents a kink at the shock position, where the characteristics cross for the first 
time. Along the dot-dashed lines R(t,x) is constant. 



The situation where > < and 9 < > is qualitatively different. In this case, the 
characteristics cross at some finite time: a shock occurs - see Fig. |l|(b). A crossing point 
of two characteristics means indeed that at this point, two different values of R (or H) are 
possible and these functions then become discontinuous. Strictly speaking, the Eqs. are 
no longer valid, and the diffusion terms left of from the analysis become important to smooth 
out this discontinuity. In Fig. |[ we plotted snapshots of the h and R profiles at different 
times, for both situations (with and without the occurrence of a shock). 

One can calculate the time t s and location x s at which the shock occurs. For that purpose, 
let us introduce the envelop of the characteristic curves x(t, v), where v is a label. The envelop 
can be represented in a parametric way as (t e (v), x e {v)). It has the property that for each 
of its points exists a characteristic, which touches it tangentially. It has then to fulfill the 
conditions x(t e (y),v) — x e [y), x v (t e (u),v) = 0. After some calculations, one can find the 
explicit expression for the envelop, 



x e (u) = v- 



10< 



1 + cr/ + A(0, v) 1 + 



1-/3 



, . . 1 A(0,j/) . / , A(0,i/) x 
07 + ln(a7) + 1 + A „ - In [ -1 - 



1-/3 



1-/3 



(11) 



(12) 
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Fig. 2. - Total height profiles h(t, x) for the two cases of Fig. [y. The bold line in (a) shows the critical 
slope, which is chosen here as 9 r = 0.2. In (b) the shock occurs at t a = 73.78, x 3 = —16.05. Insets: 
Corresponding evolution of the amount of moving grains R(t,x). For t = 75 the solutions h(t,x) and 
R(t,x) are not single valued for —15.61 < x < —15.20 since each point within the range bounded by 
the envelop is covered three times by a characteristic C_ . As orientation guide, the limiting values at 
both boundaries are connected here by straight lines. 



This envelop has two branches, separated by a kink, see Fig. |l|(b), given by v = v s = 
[cry + In (017) — 1 + f3 — ln(l — j3)]/(7#<). Whereas the upper branch is parameterized by 
—00 < v < v s , the lower one corresponds to v s < v < v c = [f3 + cry + ln(— cry//3)]/(7#<). The 
resulting shock coordinates are 



1 

i0l 



ln(cry) - ln(l -/?) + !+ ' 



1-/3 



1 



l-|)lu(1 -•/) - hl(n,) 



(13) 

The condition that (t s ,x s ) has to be located inside region II leads to the boundary between 
the classes with and without shock as mentioned above. At the shock position, the amount of 
moving grains is universal (independent of the initial value a), and given by R s = 1/(7(1 —/?)), 
while H s = 9<y s . Since typically v ~ 7c? with d the grain diameter, we have in our rescaled 
units 7 ~ 1 showing that due to R s ~ 1 non linear saturation terms can be neglected at the 
shock if /3 <~ —1. The lower branch of the envelop saturates for large t exponentially fast with 
a characteristic time l/(7#>) at = [1 +\n{—a r y/(3)]/{ / j6 < ), which is always larger than x s . 
This means that the shock stops propagating upwards. A large time expansion in the shock free 
range — t < x < x^ gives, taking the two leading terms of W, u(t,x) — —(a/8 < )exp[-f6 < t — 
{6 < / a)xe~^ 9<t ]. Thus the slope is non monotonous within this range: after increasing for 
small times it relaxes again to the initial value 8 < as H x (t,x) — 9 < exp[— (9 < /a)xe r9<t }. 

Discussion. — Let us summarize the major results of this paper, which could be explored 
experimentally. Starting from an initial profile made up of two different slopes, we find that 
shocks can occur after a finite time, depending on the value of the two slopes and the initial 
density of rolling grains. When shocks are absent, we find that the evolution surface profile is 
characterized by different velocities: the kink moves upwards with a velocity of the order of 
07 for early times, while the edge of the "active" region moves downwards at a velocity which 
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only depends on the initial slopes, and is smaller than the velocity of the grains. The final 
slope is shown to be the angle of repose; however, for finite size systems, one expects the final 
slope to be smaller by an amount which varies as 1/lnL. When a shock appears, we predict 
the time and position of this shock, as well as the density of rolling grains there, which takes 
a universal value. The shock is found to stop progressing upwards. 

Our results are in disagreement with those of MP. For the situation considered here, they 
predict that the initial profile is rigidly shifted along straight characteristics. Therefore, for 
example, the final slope would be given by H x (t,x = 0) = 9 < , which is completely different 
from our prediction of a decaying slope. The reason for this discrepancy comes from their 
implicit assumption that Ro(x) + H (x) + \n(R (x)) / j = const., which does not hold in the 
cases considered here. 

The method presented here can be extended to more general situations. For example, each 
profile Hq(x) can be approximated by a piecewise linear function. Therefore, our analysis 
can be used to obtain analytical results for more complicated situations as, e.g., bumps or 
sinusoidal shapes. Another interesting situation is the case where Ro(x) is localized in space. 
Applications of this method to the problem of ripple formation arc under way. 

Two important physical phenomena have been neglected: diffusion terms, which are ex- 
pected to be important in the presence of shocks or in the case of a localized initial Rq(x) (see 
[||), and non-linear effects, which lead to a saturation of the static/rolling grains conversion 
term. A simple way to account for the latter effect is to replace the characteristics by straight 
lines of velocity 7-Roo as soon as R = R^ . The influence of a dependence of the velocity of 
grains on their density would also be worth investigating . 

This research was partly supported by the Deutsche Forschungsgemeinschaft (DFG) under 
grant EM70/1-1. 
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